{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a target=\"_blank\" href=\"https://colab.research.google.com/github/opensim-org/opensim-core/blob/main/Bindings/Python/tutorials/Tutorial%206%20-%20Static%20Optimization.ipynb\">\n",
    "  <img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/>\n",
    "</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JLgpKeJ4NYHz"
   },
   "source": [
    "# 6. Static Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SwX6QwGU0Jwj"
   },
   "source": [
    "## 6.1. Objectives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4p_3qlB60Np3"
   },
   "source": [
    "**Introduction**\n",
    "\n",
    "Static Optimization is a method for estimating muscle activations and muscle forces that satisfy the positions, velocities, accelerations, and external forces (e.g., ground reaction forces) of a motion. The technique is called  \"static\" since calculations are performed at each time frame, without integrating the equations of motion between time steps. Static Optimization can be very fast and efficient because there is no integration, but it does ignore activation dynamics and tendon compliance. See [[1]](https://doi.org/10.1115/1.4029304) for more details regarding this and similar modeling and simulation choices and their pros and cons."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_cj6YClyLG_C"
   },
   "source": [
    "\n",
    "**Purpose**\n",
    "\n",
    "In this tutorial, you will learn how to perform a Static Optimization analysis on a single stride gait trial and assess gastrocnemius activation and force production in a simplified example model ([gait10dof18musc](https://simtk-confluence.stanford.edu:8443/display/OpenSim/Musculoskeletal+Models)).\n",
    "\n",
    "As with any analysis or simulation, the Static Optimization results' quality largely depends on the quality of the model, motion, and forces used as inputs. The model should have mass, anthropometry, and strength that represent the experimental participant; while the degrees of freedom and muscle geometry should be appropriate for the questions being asked (e.g., studying an upper body motion would require a more detailed upper extremity muscle set and skeletal geometry). The motion should contain smooth, realistic accelerations and all external forces during the motion (e.g., ground reaction forces) should be accurately measured and applied to the model. Additional forces (i.e., reserves and residuals) are often needed but should be small enough to not confound analysis. Any issues with your inputs will give you poor results or cause the tool to fail.\n",
    "\n",
    "In this tutorial, you will:\n",
    "\n",
    "*   Get an initial Static Optimization analysis running without error messages or crashes.\n",
    "*   Improve the activation and force results by using improved motion data.\n",
    "*   Learn how to reduce actuation from reserve and residual actuators.\n",
    "*   Add passive forces and observe their effect on muscle activation.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "G7_bQ5HTLJFl"
   },
   "source": [
    "**Format**\n",
    "\n",
    "Each section of the tutorial guides you through certain tools within OpenSim and asks you to answer a few questions. The questions can be answered based on information from OpenSim and basic knowledge of the human musculoskeletal system. As you complete each section of the tutorial, feel free to explore OpenSim and the lower extremity model further on your own."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "uoMuDaWxLMW1"
   },
   "source": [
    "\n",
    "**Before You Get Started**\n",
    "\n",
    "This tutorial begins in a workflow where the model is scaled, Inverse Kinematics run, and Residual Reduction Analysis performed.\n",
    "\n",
    "Before completing the tutorial, make sure you have reviewed the following reading from the User's Guide and completed the prerequisite tutorials.\n",
    "\n",
    "User Guide | Prerequisite Tutorials | Other Resources\n",
    "--------------------|------------------------|------------------\n",
    "[Static Optimization](https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53090088/Static+Optimization) | [Tutorial 1: Introduction to OpenSim](https://colab.research.google.com/github/opensim-org/opensim-core/blob/main/Bindings/Python/tutorials/Tutorial%201%20-%20Introduction%20to%20OpenSim.ipynb)          | [[1]](https://doi.org/10.1115/1.4029304)  \n",
    "[Inverse Kinematics](https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53090037/Inverse+Kinematics)  | [Tutorial 2 : Creating and Simulating a Simple Arm Model ](https://colab.research.google.com/github/opensim-org/opensim-core/blob/main/Bindings/Python/tutorials/Tutorial%202%20-%20Creating%20and%20Simulating%20a%20simple%20arm%20Model.ipynb)         | [[2]](https://doi.org/10.1109/TBME.2007.901024)  \n",
    "[Residual Reduction Algorithm](https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53089669/Residual+Reduction+Algorithm)  | [Tutorial 3: Loading and Modifying OpenSim Models](https://colab.research.google.com/github/opensim-org/opensim-core/blob/main/Bindings/Python/tutorials/Tutorial%203%20-%20Loading%20and%20Modifying%20OpenSim%20Models.ipynb)        |               \n",
    "[Inverse Dynamics](https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53090074/Inverse+Dynamics)  | [Tutorial 4: Musculoskeletal Models, Motion Files, and MuscleAnalysis](https://colab.research.google.com/github/opensim-org/opensim-core/blob/main/Bindings/Python/tutorials/Tutorial%204%20-%20Musculoskeletal%20models%2C%20Motion%20Files%20and%20MuscleAnalysis.ipynb)           |      |  \n",
    "[Scaling](https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53090000/Scaling)       | [Tutorial 5: Scaling, Inverse Kinematics, and Inverse Dynamics](https://colab.research.google.com/github/opensim-org/opensim-core/blob/main/Bindings/Python/tutorials/Tutorial%205%20-%20Scaling%2C%20Inverse%20Kinematics%2C%20and%20Inverse%20Dynamics.ipynb)           |                    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NeAzWgeVLOJL"
   },
   "source": [
    "**Resources**\n",
    "\n",
    "All of the files necessary to complete this tutorial are available in the following Google Drive folder: https://drive.google.com/drive/folders/1XFkFE018DLXxPto6XdVRcY8Fl6hP3ecA?usp=sharing\n",
    "\n",
    "While you complete the tutorial, you will find cells that automatically download the required files.\n",
    "\n",
    "To obtain a link to one of the files in Google Drive, right-click over it and click on **Get Link**.\n",
    "\n",
    "**Local Installation**\n",
    "\n",
    "You may also run this notebook using a locally installed OpenSim `conda` package. Visit [the OpenSim python scripting Confluence page](https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53085346/Scripting+in+Python) for more information on how to install OpenSim in a local Python environment. You will also need to download the resources from Google Drive to the same directory as this notebook. Then, you may skip to Section 6.3 below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ADllcxbsjTr8"
   },
   "source": [
    "## 6.2. Setup `conda` and OpenSim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mmmTuEI4HQJG"
   },
   "source": [
    "First, set up the environment by executing the following cell (See [Tutorial 1: Introduction to OpenSim](https://colab.research.google.com/github/opensim-org/opensim-core/blob/main/Bindings/Python/tutorials/Tutorial%201%20-%20Introduction%20to%20OpenSim.ipynb) for more information)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "0w9LYQXx5TO4"
   },
   "outputs": [],
   "source": [
    "!pip install -q condacolab\n",
    "import condacolab\n",
    "condacolab.install()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vt5loGbdfKb8"
   },
   "source": [
    "Don't worry if after executing the previous cell you get an error saying that your session has failed. This happens because `condacolab` needs to restart the session for the changes to take effect. Therefore, you have to execute the previous cell before executing any other cell.\n",
    "\n",
    "Now, execute the following cell to install the OpenSim `conda` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "vWf9ObZJ5TO5"
   },
   "outputs": [],
   "source": [
    "!conda install opensim-org::opensim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "56GWGq3HSLNQ"
   },
   "source": [
    "## 6.3. Study 1: Run Static Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2gJ1cEsvSUIB"
   },
   "source": [
    "The objective of this study is to perform Static Optimization over the `gait10dof18musc` model, using results from IK as motion and a set of external forces applied to it. You don't have to worry about calculating the IK results and the external forces, since you will download them in the next subsection. Once Static Optimization has been calculated, we want to visualize the activation of the gastrocnemius muscle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "337hUq4BSNHQ"
   },
   "source": [
    "**Printing OpenSim logs**\n",
    "\n",
    "To visualize OpenSim's log information, we need to redirect it to either the console, or a file. In this case, we are going to do both."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mLyWh3uuZM1W"
   },
   "source": [
    "To output OpenSim's log information to the console, we need to create a [LogSink](https://simtk.org/api_docs/opensim/api_docs/classOpenSim_1_1LogSink.html). A `LogSink` is an object that is able to report logging messages.  We want the log messages to appear as python output, so we implement the class in a way that every message is printed, and we add the sink to OpenSim's [Logger](https://simtk.org/api_docs/opensim/api_docs/classOpenSim_1_1Logger.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lu_RSBjaNXkq"
   },
   "source": [
    "To output OpenSim's log information to a file, we just need to call the function `addSinkFile()`, passing as an argument a name for the log file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kZO-sJ3sYGpG"
   },
   "outputs": [],
   "source": [
    "import opensim as osim\n",
    "\n",
    "class PythonOutputSink(osim.common.LogSink):\n",
    "    def __init__(self):\n",
    "      super().__init__()\n",
    "      self.msg = \"\"\n",
    "\n",
    "    def flushImpl(self):\n",
    "      self.msg = \"\"\n",
    "\n",
    "    def sinkImpl(self, str):\n",
    "      self.msg += str + \"\\n\";\n",
    "      print(self.msg)\n",
    "\n",
    "python_output_sink = PythonOutputSink()\n",
    "\n",
    "osim.Logger.addSink(python_output_sink)\n",
    "osim.Logger.addFileSink(\"opensim-log.log\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "e2kJm4w-YuMF"
   },
   "source": [
    "**Download the files**\n",
    "\n",
    "First, let's download the model (`gait10dof18musc_simbody.osim`), the IK results (`subject01_walk_IK.mot`) and the external forces (`subject01_walk_grf.xml` and `subject01_walk_grf.mot`) files. The files are stored in a Google Drive folder (See **Resources** in the Objectives section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "S_LJjt3YBDiZ"
   },
   "outputs": [],
   "source": [
    "!gdown \"1CB6-XNDMjNJ_Ca57AejgolT4iQTqVMGg&confirm=t\" # gait10dof18musc_simbody.osim\n",
    "!gdown \"1sPgI8VzLvp6RC9KZMgAc4fA7S4CC6jQx&confirm=t\" # subject01_walk_IK.mot\n",
    "!gdown \"1n7BN9rJv_OPW20HtFdiJr--2aqkg5Scq&confirm=t\" # subject01_walk_grf.xml\n",
    "!gdown \"1wkMVlCN7ka4Bt-x8aBVj_c5oK9rH1tW6&confirm=t\" # subject01_walk_grf.mot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bjlOR8fym2Pl"
   },
   "source": [
    "Load the model into OpenSim. To check that the file has been downloaded correctly, let's print its name. Notice that we are getting more information printed to the console because of the Python `LogSink` that we've added. You can also see the logs in the `opensim-log.log` file that has been created.\n",
    "\n",
    "Don't worry about logged messages similar to \"Couldn't find file 'sacrum.vtp'\"; VTP files are only used for visualization and not needed for this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "oGTQE8Azh6p8"
   },
   "outputs": [],
   "source": [
    "# Load the model.\n",
    "gait1018 = osim.Model('gait10dof18musc_simbody.osim')\n",
    "print(\"Name of the model:\", gait1018.getName())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GZ-Wuh6ACWpN"
   },
   "source": [
    "**Create Static Optimization Analysis Tool**\n",
    "\n",
    "\n",
    "The first step to perform a Static Optimization Analysis is to create a Static Optimization object and configure it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4yiF1t4mCRrH"
   },
   "outputs": [],
   "source": [
    "# Define a StaticOptimization object.\n",
    "so = osim.StaticOptimization()\n",
    "\n",
    "# Set start and end times for the analysis.\n",
    "start_time = 0.3\n",
    "end_time = 1.5\n",
    "so.setStartTime(start_time)\n",
    "so.setEndTime(end_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iVvH2t6sMc7t"
   },
   "source": [
    "Then, we need to create an `AnalyzeTool`. We attach two files to the `AnalyzeTool`:\n",
    "*  The IK results file (`subject01_walk_IK.mot`).\n",
    "*  The external loads file (`subject01_walk_grf.xml`).\n",
    "\n",
    "The tool will step through each time frame of the input motion. At each time, it will update the model coordinates, external loads, and calculate the muscle forces and activations needed to generate the model's accelerations based on the input motion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "YxX1dIKYEneO"
   },
   "outputs": [],
   "source": [
    "# Create analyze tool for static optimization.\n",
    "so_analyze_tool = osim.AnalyzeTool()\n",
    "so_analyze_tool.setName(\"SO\")\n",
    "\n",
    "# Set model file, motion files and external load file names.\n",
    "so_analyze_tool.setModelFilename(\"gait10dof18musc_simbody.osim\")\n",
    "so_analyze_tool.setCoordinatesFileName(\"subject01_walk_IK.mot\")\n",
    "so_analyze_tool.setExternalLoadsFileName(\"subject01_walk_grf.xml\")\n",
    "\n",
    "# Add analysis.\n",
    "so_analyze_tool.updAnalysisSet().cloneAndAppend(so)\n",
    "\n",
    "# Configure analyze tool.\n",
    "so_analyze_tool.setReplaceForceSet(False)\n",
    "so_analyze_tool.setStartTime(start_time)\n",
    "so_analyze_tool.setFinalTime(end_time)\n",
    "\n",
    "# Directory where results are stored.\n",
    "so_analyze_tool.setResultsDir(\"SO_Results\")\n",
    "\n",
    "# Print configuration of analyze tool to a xml file.\n",
    "so_analyze_tool.printToXML(\"SO_AnalyzeTool_setup.xml\")\n",
    "\n",
    "# Load configuration and run the analyses. \n",
    "so_analyze_tool = osim.AnalyzeTool(\"SO_AnalyzeTool_setup.xml\", True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6u4xoyHdQXBq"
   },
   "source": [
    "**Perform Static Optimization**\n",
    "\n",
    "With the Analysis Tool defined and configured, we can perform Static Optimization. This may take a few seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "aUbfhgubQXhu"
   },
   "outputs": [],
   "source": [
    "# Run static optimization.\n",
    "so_analyze_tool.run();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tpuoaR2SMl2l"
   },
   "source": [
    "**Plot results**\n",
    "\n",
    "With the Static Optimization execution finished, now you can extract the activation of the gastrocnemius muscle from the generated result files and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "y4CjW5RNF3RR"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_activation = osim.TimeSeriesTable(\n",
    "    \"SO_Results/SO_StaticOptimization_activation.sto\")\n",
    "\n",
    "# Get values.\n",
    "time = table_activation.getIndependentColumn()\n",
    "gastroc_activation = table_activation.getDependentColumn(\"gastroc_r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "w4J1bJGpzC3U"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Gastrocnemius Activation')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time, gastroc_activation.to_numpy(), label='Gastrocnemius (right)')\n",
    "axs.set_xlabel('Time (s)')\n",
    "axs.set_ylabel('Activation')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QGiGKzOdNIn1"
   },
   "source": [
    "**Results Explanation**\n",
    "\n",
    "What exactly is happening? If you check the output of the Static Optimization execution (when you executed run on the AnalyzeTool), there is an error message saying **Optimizer Failed...**. This is informing you that the optimizer could not find an adequate solution at the time indicated. However, it is not crashing since the tool is trying to find adequate solutions in the following times. Remember that Static Optimization is trying to calculate activations and muscle forces that satisfy the model's coordinate accelerations calculated from the specified kinematics file. If the model's actuators cannot generate sufficient forces to match these accelerations, this \"Optimizer Failed...\" error will appear. The Static Optimization tool has produced a set of results to help with troubleshooting, but these outputs are not valid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jkI_ssdCoZrM"
   },
   "source": [
    "An under-actuated model is the product of muscles not being strong enough to generate the required forces and/or model coordinates not having any actuation (muscles or other actuators). In this particular model, there are no muscles that control lumbar flexion/extension, so all the other lower limb muscles are trying to control trunk motion, without much success. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WZhmZfCeobei"
   },
   "source": [
    "One possible solution is to add extra muscles that cross the lumbar joint. However in this series of studies we are only interested in gastrocnemius activation during walking. Adding additional muscles will add extra complexity to the model that is likely not required for our analysis. Instead of lumbar muscles, we will add a torque actuator at the lumbar joint."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "uEq4VwtQodXH"
   },
   "source": [
    "We will also add reserve and residual actuators to the model. Reserve actuators can add extra actuation during portions of the gait cycle where muscles are not able to generate sufficient accelerations (e.g., during a spike in acceleration). Residual actuators are \"hand of God\" forces that account for (what should be small) discrepancies between the model, measured motions, and forces. In other words, these actuators ensure that Newton's 2$^{nd}$ law, $F=m*a$, is satisfied throughout the analysis. For more details, see _Section 3.1.3_ of [[1](https://doi.org/10.1115/1.4029304)]. Typically, these residual actuators are added to the model's body that is connected to ground (in this particular case, residual actuators will be added to the pelvis body)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-lsgJYW3Sd5W"
   },
   "source": [
    "### 6.3.1. Adding Actuators to the Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AA_lennaRyp0"
   },
   "source": [
    "You will now add a set of additional residual actuators to fix the errors obtained in the previous section. This includes a lumbar actuator, reserve actuators at the other joints, and residual actuators at the pelvis. Once you have added the actuatiors, you will rerun the analysis to check if the problem is solved."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EZ41sHvFSFZq"
   },
   "source": [
    "**Download the files**\n",
    "\n",
    "First, let's download the actuators file (`gait10dof18musc_Strong_actuators.xml`). The file is stored in a Google Drive folder (See **Resources** in the Objectives section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "9d2gSzuYFaS1"
   },
   "outputs": [],
   "source": [
    "!gdown \"1qYK7JULDucOl_GN3-74OetYJ9ckKcev6&confirm=t\" # gait10dof18musc_Strong_actuators.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iJZg28MEkPoS"
   },
   "source": [
    "**Create Static Optimization Analysis Tool**\n",
    "\n",
    "Now, let's add actuators to the Analyze Tool. You don't need to add the actuators manually since you have downloaded a file defining them (`gait10dof18musc_Strong_actuators.xml`). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AZRgXMkMof3x"
   },
   "source": [
    "In order to preserve previous results, we are going to change the name of the analysis tool and the name of the output folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Ziflsp26FAyW"
   },
   "outputs": [],
   "source": [
    "# Create analyze tool for static optimization.\n",
    "so_strong_analyze_tool = osim.AnalyzeTool()\n",
    "so_strong_analyze_tool.setName(\"SO_Strong\")\n",
    "\n",
    "# Add actuators file.\n",
    "forceSet_files = osim.ArrayStr()\n",
    "forceSet_files.append(\"gait10dof18musc_Strong_actuators.xml\")\n",
    "so_strong_analyze_tool.setForceSetFiles(forceSet_files)\n",
    "\n",
    "# Set model file, motion files and external load file names.\n",
    "so_strong_analyze_tool.setModelFilename(\"gait10dof18musc_simbody.osim\")\n",
    "so_strong_analyze_tool.setCoordinatesFileName(\"subject01_walk_IK.mot\")\n",
    "so_strong_analyze_tool.setExternalLoadsFileName(\"subject01_walk_grf.xml\")\n",
    "\n",
    "# Add analysis.\n",
    "so_strong_analyze_tool.updAnalysisSet().cloneAndAppend(so)\n",
    "\n",
    "# Configure analyze tool.\n",
    "so_strong_analyze_tool.setReplaceForceSet(False)\n",
    "so_strong_analyze_tool.setStartTime(start_time)\n",
    "so_strong_analyze_tool.setFinalTime(end_time)\n",
    "\n",
    "# Directory where results are stored.\n",
    "so_strong_analyze_tool.setResultsDir(\"SO_Strong_Results\")\n",
    "\n",
    "# Print configuration of analyze tool to a xml file.\n",
    "so_strong_analyze_tool.printToXML(\"SO_Strong_AnalyzeTool_setup.xml\")\n",
    "\n",
    "# Load configuration and run the analyses. \n",
    "so_strong_analyze_tool = osim.AnalyzeTool(\"SO_Strong_AnalyzeTool_setup.xml\", True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "O_ZXXdlfTpE_"
   },
   "source": [
    "**Perform Static Optimization**\n",
    "\n",
    "Finally, we can perform Static Optimization on the model with the actuators. Notice how this time the calculations are faster, solutions are found at each time point, and there are no errors during the execution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "OuRQXLZPTpFA"
   },
   "outputs": [],
   "source": [
    "# Run static optimization.\n",
    "so_strong_analyze_tool.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "c0wrCHCkUB-5"
   },
   "source": [
    "**Plot Results**\n",
    "\n",
    "First, plot the lumbar actuator reserves from the output file `Static_Optimization_Strong_StaticOptimization_force.sto`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "G5PXag0uGuAa"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_forces = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_Results/SO_Strong_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_force = table_forces.getIndependentColumn()\n",
    "lumbar_extension_reserve = table_forces.getDependentColumn(\"lumbar_extension_reserve\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4Tu63Acz2wBQ"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Lumbar Reserve')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time_force, lumbar_extension_reserve.to_numpy(), label='Lumbar Reserve')\n",
    "axs.set_xlabel('Time')\n",
    "axs.set_ylabel('Moment (Nm)')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fENfqQHlUSNo"
   },
   "source": [
    "Now, let's compare the gastrocnemius muscle activation curves from the two analyses (with and without additional actuators)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ND4iTk6OJn6P"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_activation = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_Results/SO_Strong_StaticOptimization_activation.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_actuation = table_activation.getIndependentColumn()\n",
    "gastroc_activation_actuation = table_activation.getDependentColumn(\"gastroc_r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "RA4njbtC27cj"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Gastrocnemius Activation')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time, gastroc_activation.to_numpy(), label='Without Actuation')\n",
    "axs.plot(time_actuation, gastroc_activation_actuation.to_numpy(), label='With Actuation')\n",
    "axs.set_xlabel('Time (s)')\n",
    "axs.set_ylabel('Activation')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IJlwI5csUbnj"
   },
   "source": [
    "The previous figure shows that the gastrocnemius activation was significantly reduced. The reduction is due to the additional torques generated by the actuators added to the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AHQ0fSvqUhWj"
   },
   "source": [
    "Plot the moments produced by the right ankle reserve actuator from the output file `Static_Optimization_Strong_StaticOptimization_force.sto`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "51pJOOBCJ402"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_force = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_Results/SO_Strong_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_force = table_force.getIndependentColumn()\n",
    "ankle_angle_r_reserve_force = table_force.getDependentColumn(\"ankle_angle_r_reserve\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fvVWUXMu3Jge"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Ankle Reserve Moment')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time_actuation, ankle_angle_r_reserve_force.to_numpy(), label='Ankle Reserve')\n",
    "axs.set_xlabel('Time')\n",
    "axs.set_ylabel('Moment (Nm)')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9bKfj79eUqhu"
   },
   "source": [
    "Finally, plot the forces `Fx`, `Fy`, and moment `Mz` residual actuator forces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "mqilA22pKe_n"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_force = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_Results/SO_Strong_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_force = table_force.getIndependentColumn()\n",
    "fx = table_force.getDependentColumn(\"FX\")\n",
    "fy = table_force.getDependentColumn(\"FY\")\n",
    "mz = table_force.getDependentColumn(\"MZ\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "J65OKGUI3UnJ"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Residual Forces and Moment')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time_force, fx.to_numpy(), label='FX')\n",
    "axs.plot(time_force, fy.to_numpy(), label='FY')\n",
    "axs.plot(time_force, mz.to_numpy(), label='MZ')\n",
    "axs.set_xlabel('Time')\n",
    "axs.set_ylabel('Forces (N; FX, FY) & Moment (Nm; MZ)')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ppWz47WEUxt6"
   },
   "source": [
    "**Results Explanation:**\n",
    "\n",
    "We now have a working Static Optimization analysis. A functioning analysis does not mean that the results are valid; rather, we have a place to start iterating from to improve the output. In the following studies, we will gradually improve the quality of the inputs to get a better result. Review the questions below to better understand the problems with the current set of results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "R5LoTvT-ipdJ"
   },
   "source": [
    "**Questions:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "okZLnDQ-mTuG"
   },
   "source": [
    "*   The activation and force results are very noisy; why do you think this is? We will learn how to fix this in Study 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QVA_u_pTmVlz"
   },
   "source": [
    "*   The Mz contribution is between ±50 N-m. What is the contribution of the residual actuator force for Fx, and Fy? How large are these residual actuator forces in comparison to the peak ground reaction forces measured doing the motion? What improvements could you make to any of the inputs (model, kinematics, external forces) to reduce the size of the residual forces? Study 2 will address these residual results. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fpknvfWjmXPi"
   },
   "source": [
    "*   The peak ankle extension moment from Inverse Dynamics was 120 N-m. How does this peak ankle moment compare to the peak moment generated by the ankle reserve actuator (e.g., as a percentage)? Do you think the reserve actuator peak should be this large? How might we address this issue? We will learn how to fix this in Study 3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "J6_F0oT8WlVA"
   },
   "source": [
    "## 6.4. Study 2: Improving the input kinematics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pgmwSu383oCu"
   },
   "source": [
    "The aim of Study 2 is to assess how improving the input kinematics affects muscle activations and forces calculated by Static Optimization. Static Optimization is trying to compute muscle forces that generate the accelerations seen in the input kinematics. If these accelerations are noisy, there will be noise in the muscle activations and forces. In Study 1, we saw noisy activations and forces that are most likely the result of noisy motion data. We will investigate two methods for improving the motion data:\n",
    "*  Filtering the input kinematics.\n",
    "*  Using the Residual Reduction Algorithm (RRA) to generate a dynamically consistent motion with consistent accelerations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VG0ylYq_j8Ok"
   },
   "source": [
    "Filtering the kinematics will spline-fit the motion data and filter the coordinate positions, reducing noise in the accelerations. If you use IK results as inputs to Static Optimization, you should always filter them, either externally (e.g., with Matlab or Python) or by using the OpenSim filters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oWf-l1v0j-VF"
   },
   "source": [
    "Since RRA uses a forward dynamics simulation (i.e., the model's equations of motion are integrated forward in time), the output kinematics will be consistent accelerations that are great to use in Static Optimization. If you are using RRA results, you don't need to do any additional filtering. Using RRA also has the added benefit of generating a more \"dynamically consistent\" model and set of kinematics and forces. We are not going to go in depth into RRA in this particular tutorial so if you are unfamiliar with RRA please refer to the [User's Guide](https://simtk-confluence.stanford.edu:8443/display/OpenSim/Residual+Reduction+Algorithm)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UDysgMhAk9f6"
   },
   "source": [
    "### 6.4.1. Filter the motion using a low pass filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7IZDMDJLkp6N"
   },
   "source": [
    "**Create Static Optimization Analysis Tool**\n",
    "\n",
    "Let's create a Static Optimization Tool, but this time we have to append the reserve and residal actuators from the beginning, filter the motion using 6Hz cut-off, and change the output directory name, so results are not overwritten."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "b54SwBByXD-7"
   },
   "outputs": [],
   "source": [
    "# Create analyze tool for static optimization.\n",
    "so_strong_filter_analyze_tool = osim.AnalyzeTool()\n",
    "so_strong_filter_analyze_tool.setName(\"SO_Strong_Filter\")\n",
    "\n",
    "# Add actuators file.\n",
    "forceSet_files = osim.ArrayStr()\n",
    "forceSet_files.append(\"gait10dof18musc_Strong_actuators.xml\")\n",
    "so_strong_filter_analyze_tool.setForceSetFiles(forceSet_files)\n",
    "\n",
    "# Set model file, motion files and external load file names.\n",
    "so_strong_filter_analyze_tool.setModelFilename(\"gait10dof18musc_simbody.osim\")\n",
    "so_strong_filter_analyze_tool.setCoordinatesFileName(\"subject01_walk_IK.mot\")\n",
    "so_strong_filter_analyze_tool.setExternalLoadsFileName(\"subject01_walk_grf.xml\")\n",
    "\n",
    "# Change low pass frequency to filter the model.\n",
    "so_strong_filter_analyze_tool.setLowpassCutoffFrequency(6)\n",
    "\n",
    "# Add analysis.\n",
    "so_strong_filter_analyze_tool.updAnalysisSet().cloneAndAppend(so)\n",
    "\n",
    "# Configure analyze tool.\n",
    "so_strong_filter_analyze_tool.setReplaceForceSet(False)\n",
    "so_strong_filter_analyze_tool.setStartTime(start_time)\n",
    "so_strong_filter_analyze_tool.setFinalTime(end_time)\n",
    "\n",
    "# Directory where results are stored.\n",
    "so_strong_filter_analyze_tool.setResultsDir(\"SO_Strong_Filter_Results\")\n",
    "\n",
    "# Print configuration of analyze tool to a xml file.\n",
    "so_strong_filter_analyze_tool.printToXML(\"SO_Strong_Filter_AnalyzeTool_setup.xml\")\n",
    "\n",
    "# Load configuration and run the analyses. \n",
    "so_strong_filter_analyze_tool = osim.AnalyzeTool(\"SO_Strong_Filter_AnalyzeTool_setup.xml\", True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZU6sest-lTUH"
   },
   "source": [
    "**Perform Static Optimization**\n",
    "\n",
    "We can perform Static Optimization on the model with the filtered motion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DMyIhe-0XQml"
   },
   "outputs": [],
   "source": [
    "# Run static optimization.\n",
    "so_strong_filter_analyze_tool.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xZBT-ScKlesI"
   },
   "source": [
    "### 6.4.2. Smoothing the motion using RRA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ikki8Olol3J8"
   },
   "source": [
    "**Download the files**\n",
    "\n",
    "First, let's download the RRA motion file (`subject_adjusted_Kinematics_q.sto`). The file is stored in a Google Drive folder (See **Resources** in the Objectives section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "H6eAcbj6jxkV"
   },
   "outputs": [],
   "source": [
    "!gdown \"1XpgUlSGV6T7SKRnCpcIomg5aZ4INreI2&confirm=t\" # subject_adjusted_Kinematics_q.sto"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_xBCmYBZmK_H"
   },
   "source": [
    "**Create a Static Optimization Analysis Tool**\n",
    "\n",
    "Let's create a Static Optimization Tool, but this time we have to use the motion file kinematics generated by RRA (`subject_adjusted_Kinematics_q.sto`). Remember we are using RRA results so the kinematics are smoother. This time, do not apply a low pass filter, and change the step interval to 10. Again, change the output folder for the file so we don't overwrite previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "i___D9dylXCL"
   },
   "outputs": [],
   "source": [
    "# Define a StaticOptimization object.\n",
    "so_rra = osim.StaticOptimization()\n",
    "\n",
    "# Set start and end times for the analysis.\n",
    "start_time = 0.3\n",
    "end_time = 1.5\n",
    "so_rra.setStartTime(start_time)\n",
    "so_rra.setEndTime(end_time)\n",
    "\n",
    "# Step interval.\n",
    "so_rra.setStepInterval(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DTyCVAUHi8L3"
   },
   "outputs": [],
   "source": [
    "# Create analyze tool for static optimization.\n",
    "so_strong_RRA_analyze_tool = osim.AnalyzeTool()\n",
    "so_strong_RRA_analyze_tool.setName(\"SO_Strong_RRA\")\n",
    "\n",
    "# Add actuators file.\n",
    "forceSet_files = osim.ArrayStr()\n",
    "forceSet_files.append(\"gait10dof18musc_Strong_actuators.xml\")\n",
    "so_strong_RRA_analyze_tool.setForceSetFiles(forceSet_files)\n",
    "\n",
    "# Set model file, motion files and external load file names.\n",
    "so_strong_RRA_analyze_tool.setModelFilename(\"gait10dof18musc_simbody.osim\")\n",
    "so_strong_RRA_analyze_tool.setCoordinatesFileName(\"subject_adjusted_Kinematics_q.sto\")\n",
    "so_strong_RRA_analyze_tool.setExternalLoadsFileName(\"subject01_walk_grf.xml\")\n",
    "\n",
    "# Add analysis.\n",
    "so_strong_RRA_analyze_tool.updAnalysisSet().cloneAndAppend(so_rra)\n",
    "\n",
    "# Configure analyze tool.\n",
    "so_strong_RRA_analyze_tool.setReplaceForceSet(False)\n",
    "so_strong_RRA_analyze_tool.setStartTime(start_time)\n",
    "so_strong_RRA_analyze_tool.setFinalTime(end_time)\n",
    "\n",
    "# Directory where results are stored.\n",
    "so_strong_RRA_analyze_tool.setResultsDir(\"SO_Strong_RRA_Results\")\n",
    "\n",
    "# Print configuration of analyze tool to a xml file.\n",
    "so_strong_RRA_analyze_tool.printToXML(\"SO_Strong_RRA_AnalyzeTool_setup.xml\")\n",
    "\n",
    "# Load configuration and run the analyses. \n",
    "so_strong_RRA_analyze_tool = osim.AnalyzeTool(\"SO_Strong_RRA_AnalyzeTool_setup.xml\", True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dKbZ-zO2myqS"
   },
   "source": [
    "**Perform Static Optimization**\n",
    "\n",
    "We can now perform Static Optimization on the model with the smoothened motion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "J1qSYStBjn1_"
   },
   "outputs": [],
   "source": [
    "# Run static optimization.\n",
    "so_strong_RRA_analyze_tool.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Gd5hHheSm5zn"
   },
   "source": [
    "### 6.4.3. Comparing results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "uj4UPhhrlvOg"
   },
   "source": [
    "You should now have four analysis results in separate folders\n",
    "\n",
    "*  **1. SO_Results**: Static Optimization, no extra actuation (residuals and reserves), unfiltered IK motion.\n",
    "*  **2. SO_Strong_Results**: Static Optimization, with extra actuation, unfiltered IK motion.\n",
    "*  **3. SO_Strong_Filter_Results**: Static Optimization, with extra actuation, filtered IK motion.\n",
    "*  **4. SO_Strong_RRA_Results**: Static Optimization, with extra actuation, RRA motion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hOZztbjRmHJs"
   },
   "source": [
    "**Plot Results**\n",
    "\n",
    "Generate figures that compare the gastrocnemius activations for the three different kinematics inputs (**2. unfiltered**, **3. filtered**, and **4. RRA**) where reserve and residual actuators were included (i.e., skip case 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6VZffDUjmMj0"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_activation_SO_Strong = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_Results/SO_Strong_StaticOptimization_activation.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Strong = table_activation_SO_Strong.getIndependentColumn()\n",
    "gastrocnemius_muscle_activation_strong = table_activation_SO_Strong.getDependentColumn(\"gastroc_r\")\n",
    "\n",
    "# Get results from file.\n",
    "table_activation_SO_Strong_Filter = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_Filter_Results/SO_Strong_Filter_StaticOptimization_activation.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Strong_Filter = table_activation_SO_Strong_Filter.getIndependentColumn()\n",
    "gastrocnemius_muscle_activation_strong_filtered = table_activation_SO_Strong_Filter.getDependentColumn(\"gastroc_r\")\n",
    "\n",
    "# Get results from file.\n",
    "table_activation_SO_Strong_RRA = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_RRA_Results/SO_Strong_RRA_StaticOptimization_activation.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Strong_RRA = table_activation_SO_Strong_RRA.getIndependentColumn()\n",
    "gastrocnemius_muscle_activation_strong_rra = table_activation_SO_Strong_RRA.getDependentColumn(\"gastroc_r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "zpbV3yWj3_kH"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Gastrocnemius Activation')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time_SO_Strong, gastrocnemius_muscle_activation_strong.to_numpy(), label='Unfiltered')\n",
    "axs.plot(time_SO_Strong_Filter, gastrocnemius_muscle_activation_strong_filtered.to_numpy(), label='Filtered')\n",
    "axs.plot(time_SO_Strong_RRA, gastrocnemius_muscle_activation_strong_rra.to_numpy(), label='RRA')\n",
    "axs.set_xlabel('Time (s)')\n",
    "axs.set_ylabel('Activation')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xUrKLshWnzsp"
   },
   "source": [
    "Plot the `Fx` residual for the three kinematics inputs (**2. unfiltered**, **3. filtered**, and **4. RRA**)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NINOvRDwn0hD"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_forces_SO_Strong = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_Results/SO_Strong_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Strong = table_forces_SO_Strong.getIndependentColumn()\n",
    "fx_strong = table_forces_SO_Strong.getDependentColumn(\"FX\")\n",
    "\n",
    "# Get results from file.\n",
    "table_forces_SO_Strong_Filter = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_Filter_Results/SO_Strong_Filter_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Strong_Filter = table_forces_SO_Strong_Filter.getIndependentColumn()\n",
    "fx_strong_filtered = table_forces_SO_Strong_Filter.getDependentColumn(\"FX\")\n",
    "\n",
    "# Get results from file.\n",
    "table_forces_SO_Strong_RRA = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_RRA_Results/SO_Strong_RRA_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Strong_RRA = table_forces_SO_Strong_RRA.getIndependentColumn()\n",
    "fx_strong_rra = table_forces_SO_Strong_RRA.getDependentColumn(\"FX\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "wJ8l3LeY4TgY"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Fx Residual Force')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time_SO_Strong, fx_strong.to_numpy(), label='Unfiltered')\n",
    "axs.plot(time_SO_Strong_Filter, fx_strong_filtered.to_numpy(), label='Filtered')\n",
    "axs.plot(time_SO_Strong_RRA, fx_strong_rra.to_numpy(), label='RRA')\n",
    "axs.set_xlabel('Time (s)')\n",
    "axs.set_ylabel('Force (N)')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "n_EcxYFgoQeV"
   },
   "source": [
    "**Questions**:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wy2agd2_npKU"
   },
   "source": [
    "*  Describe the changes in muscle activation for the three different motion inputs. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3R2D3eWEnjxv"
   },
   "source": [
    "*  Why did the residual force reduce from the Unfiltered to Filtered condition? Why did the residual force reduce further when using the model and motion from RRA?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Rc57Z-AdnlJZ"
   },
   "source": [
    "*  Which kinematic input should be used for the Static Optimization analysis? Explain your reason(s)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bJ7pb8YEnnJZ"
   },
   "source": [
    "*  Optional: Explore some of the other forces and activation results from the analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-4q6_ZTFofwj"
   },
   "source": [
    "## 6.5. Study 3: Tuning the reserve and residual actuators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NXZyJcBtpDdU"
   },
   "source": [
    "In Study 2, we added actuators to assist each coordinate of the model. However, we saw that these reserve actuators contributed a significant amount of the total joint torque that we would prefer the muscles to provide. In this study, you will reduce the strength of the reserve actuators and observe the resulting changes in reserve actuator torques and gastrocnemius muscle activation and forces. \n",
    "\n",
    "The current actuators file has optimal forces of $100 N$ for each force actuator (i.e., residuals $FX$ and $FY$) and optimal forces of $100 Nm$ for each torque actuator (i.e., residual $MZ$, and reserve coordinate actuators). When trying to minimize the sum of activations squared, the optimizer can use small activation values to these force and torque actuators to generate large forces or torques. For example, with an optimal torque of $100 Nm$, an activation of $0.5$ will generate a torque of $50 Nm$. If, instead, the optimal force were $1 Nm$, the optimizer must now use very high activations to produce any meaningful torques. In this latter case (optimal torque = 1), an activation of $0.5$ would only generate $0.5 Nm$, and an activation of $50$ is necessary to generate $50 Nm$ of torque. Thus, when the optimal force is very low ($~1 N$ or $1 Nm$) the optimizer can still choose to use the reserve actuators when needed, but the cost is so high relative to the muscles that the optimizer will preferentially use muscles to generate joint torques over the reserve actuators."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "n5P-y18ssZQ3"
   },
   "source": [
    "**Edit Actuators File**\n",
    "\n",
    "Edit the actuators file and change each actuator's optimal force from 100 to 1. You can do this programatically by loading the XML file and modifying the values of every element with tag `optimal_force`. The following function does this for you and saves the modified file as `gait10dof18musc_Actuators_Normal.xml`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-jyIrh2IvERk"
   },
   "outputs": [],
   "source": [
    "from xml.dom import minidom\n",
    "\n",
    "# Parse the xml file.\n",
    "file = minidom.parse('gait10dof18musc_Strong_actuators.xml')\n",
    "\n",
    "# Get all elements with tag name 'optimal_force'.\n",
    "optimal_forces = file.getElementsByTagName(\"optimal_force\")\n",
    "\n",
    "# Set all optimal forces to 1.\n",
    "for optimal_force in optimal_forces:\n",
    "  optimal_force.firstChild.nodeValue = 1\n",
    "\n",
    "xml_str = file.toprettyxml(indent =\"\\t\", newl='\\r') \n",
    "\n",
    "with open(\"gait10dof18musc_Actuators_Normal.xml\", \"w\") as f:\n",
    "    f.write(xml_str) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Sa8dYkuRsy5M"
   },
   "source": [
    "**Create Static Optimization Analysis Tool**\n",
    "\n",
    "Perform a new Static Optimization Analysis, but this time use the new Additional Forces Set file (`gait10dof18musc_Actuators_Normal.xml`). Save the results in a different folder so previous results are not overwritten."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "wETSKSk9xIte"
   },
   "outputs": [],
   "source": [
    "# Create analyze tool for static optimization.\n",
    "so_normal_RRA_analyze_tool = osim.AnalyzeTool()\n",
    "so_normal_RRA_analyze_tool.setName(\"SO_Normal_RRA\")\n",
    "\n",
    "# Add actuators file.\n",
    "forceSet_files = osim.ArrayStr()\n",
    "forceSet_files.append(\"gait10dof18musc_Actuators_Normal.xml\")\n",
    "so_normal_RRA_analyze_tool.setForceSetFiles(forceSet_files)\n",
    "\n",
    "# Set model file, motion files and external load file names.\n",
    "so_normal_RRA_analyze_tool.setModelFilename(\"gait10dof18musc_simbody.osim\")\n",
    "so_normal_RRA_analyze_tool.setCoordinatesFileName(\"subject_adjusted_Kinematics_q.sto\")\n",
    "so_normal_RRA_analyze_tool.setExternalLoadsFileName(\"subject01_walk_grf.xml\")\n",
    "\n",
    "# Add analysis.\n",
    "so_normal_RRA_analyze_tool.updAnalysisSet().cloneAndAppend(so_rra)\n",
    "\n",
    "# Configure analyze tool.\n",
    "so_normal_RRA_analyze_tool.setReplaceForceSet(False)\n",
    "so_normal_RRA_analyze_tool.setStartTime(start_time)\n",
    "so_normal_RRA_analyze_tool.setFinalTime(end_time)\n",
    "\n",
    "# Directory where results are stored.\n",
    "so_normal_RRA_analyze_tool.setResultsDir(\"SO_Normal_RRA_Results\")\n",
    "\n",
    "# Print configuration of analyze tool to a xml file.\n",
    "so_normal_RRA_analyze_tool.printToXML(\"SO_Normal_RRA_AnalyzeTool_setup.xml\")\n",
    "\n",
    "# Load configuration and run the analyses. \n",
    "so_normal_RRA_analyze_tool = osim.AnalyzeTool(\"SO_Normal_RRA_AnalyzeTool_setup.xml\", True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ViLGeIUsuLcQ"
   },
   "source": [
    "**Perform Static Optimization**\n",
    "\n",
    "We can now perform Static Optimization with the residual and reserve actuators each with an optimal force of 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "esvOnBh4xsGN"
   },
   "outputs": [],
   "source": [
    "# Run static optimization.\n",
    "so_normal_RRA_analyze_tool.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xm4GhhLVuGos"
   },
   "source": [
    "**Plot Results**\n",
    "\n",
    "Once the analysis has finished, plot the gastrocnemius activation when using the strong actuators (optimal force of 100) vs the normal actuators (optimal force of 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1J4KRZnux6Xp"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_activation_SO_Strong_RRA = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_RRA_Results/SO_Strong_RRA_StaticOptimization_activation.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Strong_RRA = table_activation_SO_Strong_RRA.getIndependentColumn()\n",
    "gastrocnemius_muscle_activation_strong_rra = table_activation_SO_Strong_RRA.getDependentColumn(\"gastroc_r\")\n",
    "\n",
    "# Get results from file.\n",
    "table_activation_SO_Normal_RRA = osim.TimeSeriesTable(\n",
    "    \"SO_Normal_RRA_Results/SO_Normal_RRA_StaticOptimization_activation.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Normal_RRA = table_activation_SO_Normal_RRA.getIndependentColumn()\n",
    "gastrocnemius_muscle_activation_normal_rra = table_activation_SO_Normal_RRA.getDependentColumn(\"gastroc_r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DTDNDHKwEuoD"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Gastrocnemius Activation')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time_SO_Normal_RRA, gastrocnemius_muscle_activation_normal_rra.to_numpy(), label='Normal')\n",
    "axs.plot(time_SO_Strong_RRA, gastrocnemius_muscle_activation_strong_rra.to_numpy(), label='Strong')\n",
    "axs.set_xlabel('Time (s)')\n",
    "axs.set_ylabel('Activation')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wsqvkJZtySfS"
   },
   "source": [
    "Plot and compare the lumbar actuator moment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "OXnYt2D0yTLh"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_forces_SO_Strong_RRA = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_RRA_Results/SO_Strong_RRA_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Strong_RRA = table_forces_SO_Strong_RRA.getIndependentColumn()\n",
    "lumbar_actuator_moment_strong_rra = table_forces_SO_Strong_RRA.getDependentColumn(\"lumbar_extension_reserve\")\n",
    "\n",
    "# Get results from file.\n",
    "table_forces_SO_Normal_RRA = osim.TimeSeriesTable(\n",
    "    \"SO_Normal_RRA_Results/SO_Normal_RRA_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Normal_RRA = table_forces_SO_Normal_RRA.getIndependentColumn()\n",
    "lumbar_actuator_moment_normal_rra = table_forces_SO_Normal_RRA.getDependentColumn(\"lumbar_extension_reserve\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "h4vIMTvxHTw7"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Lumbar Actuator Moment')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time_SO_Normal_RRA, lumbar_actuator_moment_normal_rra.to_numpy(), label='Normal')\n",
    "axs.plot(time_SO_Strong_RRA, lumbar_actuator_moment_strong_rra.to_numpy(), label='Strong')\n",
    "axs.set_xlabel('Time')\n",
    "axs.set_ylabel('Moment (Nm)')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sBicC1_Hu90Y"
   },
   "source": [
    "Plot and compare the ankle reserve actuator moment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "r4SKFl0H3O8V"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_forces_SO_Strong_RRA = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_RRA_Results/SO_Strong_RRA_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Strong_RRA = table_forces_SO_Strong_RRA.getIndependentColumn()\n",
    "ankle_r_actuator_moment_strong_rra = table_forces_SO_Strong_RRA.getDependentColumn(\"ankle_angle_r_reserve\")\n",
    "\n",
    "# Get results from file.\n",
    "table_forces_SO_Normal_RRA = osim.TimeSeriesTable(\n",
    "    \"SO_Normal_RRA_Results/SO_Normal_RRA_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Normal_RRA = table_forces_SO_Normal_RRA.getIndependentColumn()\n",
    "ankle_r_actuator_moment_normal_rra = table_forces_SO_Normal_RRA.getDependentColumn(\"ankle_angle_r_reserve\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "L5wO-niVHnQL"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Ankle Actuator Moment')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time_SO_Normal_RRA, ankle_r_actuator_moment_normal_rra.to_numpy(), label='Normal')\n",
    "axs.plot(time_SO_Strong_RRA, ankle_r_actuator_moment_strong_rra.to_numpy(), label='Strong')\n",
    "axs.set_xlabel('Time')\n",
    "axs.set_ylabel('Moment (Nm)')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kIg-uimyvGb6"
   },
   "source": [
    "Finally, plot and compare the `Fx` residual force when using the strong actuators (optimal force at 100) and normal actuators (optimal force at 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "TrSzwaOL3nJN"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_forces_SO_Strong_RRA = osim.TimeSeriesTable(\n",
    "    \"SO_Strong_RRA_Results/SO_Strong_RRA_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Strong_RRA = table_forces_SO_Strong_RRA.getIndependentColumn()\n",
    "fx_strong_rra = table_forces_SO_Strong_RRA.getDependentColumn(\"FX\")\n",
    "\n",
    "# Get results from file.\n",
    "table_forces_SO_Normal_RRA = osim.TimeSeriesTable(\n",
    "    \"SO_Normal_RRA_Results/SO_Normal_RRA_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Normal_RRA = table_forces_SO_Normal_RRA.getIndependentColumn()\n",
    "fx_normal_rra = table_forces_SO_Normal_RRA.getDependentColumn(\"FX\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "69CuU10KHxEy"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Fx Residual Force')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time_SO_Normal_RRA, fx_normal_rra.to_numpy(), label='Normal')\n",
    "axs.plot(time_SO_Strong_RRA, fx_strong_rra.to_numpy(), label='Strong')\n",
    "axs.set_xlabel('Time')\n",
    "axs.set_ylabel('Force (N)')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vDOFJWSu33fZ"
   },
   "source": [
    "**Questions**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jTTT7T3Qn5ib"
   },
   "source": [
    "*  How does the gastrocnemius muscle force change between the analysis using strong actuators (optimal force at 100) and the normal (optimal force at 1) actuators? What about the moment generated by the ankle reserve actuators? If there are differences, explain why."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "H5SWd_aqn7SP"
   },
   "source": [
    "*  How does the lumbar actuator reserve moment change between using the strong vs. normal optimal force? Why might this be the case? (Hint: observe whether or not there are muscles that control the lumbar extension degree of freedom.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UuHlogCyn8t0"
   },
   "source": [
    "*  What is the peak `Fx` residual force? Why does the value stay the same between the strong and normal actuator case? (Hint: refer to the section earlier in the tutorial where we discuss why we need residual forces and moments. Are there other actuators in the model that could generate these forces?)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Zw-iFZjs4N3E"
   },
   "source": [
    "## 6.6. Study 4: Static Optimization with a passive element \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dFVokuxe4RG1"
   },
   "source": [
    "In this Study, you will add an ankle spring to the model and observe the effect on gastrocnemius activation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Iqz8EkHQvZzq"
   },
   "source": [
    "**Add a Spring to the Model**\n",
    "\n",
    "The first step is to add a coordinate spring to the model. The following cell adds the spring for you. A new model named `gait10dof18musc_spring.osim` will be generated. Explore the different properties of the spring (e.g., stiffness, limits, damping, etc.)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DuCHCy594v54"
   },
   "outputs": [],
   "source": [
    "# Get a handle to the current model and create a new copy \n",
    "ankleSpringModel = osim.Model('gait10dof18musc_simbody.osim')\n",
    "ankleSpringModel.setName(ankleSpringModel.getName()+'_ankle_spring')\n",
    "\n",
    "# Create the spring we'll add to the model (a CoordinateLimitForce in OpenSim)\n",
    "ankleSpring = osim.CoordinateLimitForce()\n",
    "\n",
    "# Set the coordinate for the spring\n",
    "ankleSpring.set_coordinate('ankle_angle_r')\n",
    "ankleSpring.setName('AnkleLimitSpringDamper')\n",
    "\n",
    "# Add the spring to the model\n",
    "ankleSpringModel.addForce(ankleSpring)\n",
    "\n",
    "# Set the spring's properties\n",
    "ankleSpring.setUpperStiffness(30.0)\n",
    "ankleSpring.setUpperLimit(5.0)\n",
    "ankleSpring.setLowerStiffness(1.0)\n",
    "ankleSpring.setLowerLimit(-90.0)\n",
    "ankleSpring.setDamping(0.01)\n",
    "ankleSpring.setTransition(2.0)\n",
    "\n",
    "# Save the model to file\n",
    "fullPathName = ankleSpringModel.getInputFileName()\n",
    "newName = fullPathName.replace('.osim', '_spring.osim')\n",
    "ankleSpringModel.printToXML(newName)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sjZ6ynSh5sbg"
   },
   "source": [
    "**Create Static Optimization Analysis Tool**\n",
    "\n",
    "Now that you have a new model with an ankle spring, you will re-run the Static Optimization analysis using the same setup from Study 3 (with normal reserve actuators). Rename the output folder to `\\SO_Normal_RRA_Spring_Results`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "wL5Qr7NO57Jb"
   },
   "outputs": [],
   "source": [
    "# Create analyze tool for static optimization.\n",
    "so_normal_RRA_spring_analyze_tool = osim.AnalyzeTool()\n",
    "so_normal_RRA_spring_analyze_tool.setName(\"SO_Normal_RRA_Spring\")\n",
    "\n",
    "# Add actuators file.\n",
    "forceSet_files = osim.ArrayStr()\n",
    "forceSet_files.append(\"gait10dof18musc_Actuators_Normal.xml\")\n",
    "so_normal_RRA_spring_analyze_tool.setForceSetFiles(forceSet_files)\n",
    "\n",
    "# Set model file, motion files and external load file names.\n",
    "so_normal_RRA_spring_analyze_tool.setModelFilename(\"gait10dof18musc_simbody_spring.osim\")\n",
    "so_normal_RRA_spring_analyze_tool.setCoordinatesFileName(\"subject_adjusted_Kinematics_q.sto\")\n",
    "so_normal_RRA_spring_analyze_tool.setExternalLoadsFileName(\"subject01_walk_grf.xml\")\n",
    "\n",
    "# Add analysis.\n",
    "so_normal_RRA_spring_analyze_tool.updAnalysisSet().cloneAndAppend(so_rra)\n",
    "\n",
    "# Configure analyze tool.\n",
    "so_normal_RRA_spring_analyze_tool.setReplaceForceSet(False)\n",
    "so_normal_RRA_spring_analyze_tool.setStartTime(start_time)\n",
    "so_normal_RRA_spring_analyze_tool.setFinalTime(end_time)\n",
    "\n",
    "# Directory where results are stored.\n",
    "so_normal_RRA_spring_analyze_tool.setResultsDir(\"SO_Normal_RRA_Spring_Results\")\n",
    "\n",
    "# Print configuration of analyze tool to a xml file.\n",
    "so_normal_RRA_spring_analyze_tool.printToXML(\"SO_Normal_RRA_Spring_AnalyzeTool_setup.xml\")\n",
    "\n",
    "# Load configuration and run the analyses. \n",
    "so_normal_RRA_spring_analyze_tool = osim.AnalyzeTool(\"SO_Normal_RRA_Spring_AnalyzeTool_setup.xml\", True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "80i0nMGmwHki"
   },
   "source": [
    "**Perform Static Optimization**\n",
    "\n",
    "We can now perform Static Optimization on the model with the ankle spring."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "r_rlitXC57Jc"
   },
   "outputs": [],
   "source": [
    "# Run static optimization.\n",
    "so_normal_RRA_spring_analyze_tool.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "elsHoKnQ6rmg"
   },
   "source": [
    "**Plot Results**\n",
    "\n",
    "Plot the AnkleLimitSpringDamper moment from the file `SO_Normal_RRA_Spring_Results/SO_Normal_RRA_Spring_StaticOptimization_force.sto`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "lroKJsWX6wYp"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_forces_SO_Normal_RRA_Spring = osim.TimeSeriesTable(\n",
    "    \"SO_Normal_RRA_Spring_Results/SO_Normal_RRA_Spring_StaticOptimization_force.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Strong_RRA = table_forces_SO_Normal_RRA_Spring.getIndependentColumn()\n",
    "ankle_limit_spring_damper = table_forces_SO_Normal_RRA_Spring.getDependentColumn(\"AnkleLimitSpringDamper\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "STLaD20tH93L"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Ankle Coordinate Spring')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time_SO_Strong_RRA, ankle_limit_spring_damper.to_numpy(), label='Ankle Spring')\n",
    "axs.set_xlabel('Time (s)')\n",
    "axs.set_ylabel('Ankle Limit Spring Damper')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RqkdppgY7ygX"
   },
   "source": [
    "Plot the Gastrocnemius activation from the analyzes with and without the Spring (the results should be in the folder `SO_Normal_RRA_Results` from Study 4 and the latest folder, `SO_Normal_RRA_Spring_Results`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "PGI07Ea572rf"
   },
   "outputs": [],
   "source": [
    "# Get results from file.\n",
    "table_activation_SO_Normal_RRA_Spring = osim.TimeSeriesTable(\n",
    "    \"SO_Normal_RRA_Spring_Results/SO_Normal_RRA_Spring_StaticOptimization_activation.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Normal_RRA_Spring = table_activation_SO_Normal_RRA_Spring.getIndependentColumn()\n",
    "gastrocnemius_muscle_activation_normal_rra_spring = table_activation_SO_Normal_RRA_Spring.getDependentColumn(\"gastroc_r\")\n",
    "\n",
    "# Get results from file.\n",
    "table_activation_SO_Normal_RRA = osim.TimeSeriesTable(\n",
    "    \"SO_Normal_RRA_Results/SO_Normal_RRA_StaticOptimization_activation.sto\")\n",
    "\n",
    "# Get values.\n",
    "time_SO_Normal_RRA = table_activation_SO_Normal_RRA.getIndependentColumn()\n",
    "gastrocnemius_muscle_activation_normal_rra = table_activation_SO_Normal_RRA.getDependentColumn(\"gastroc_r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "eCXPMru9IKn5"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Gastrocnemius Activation')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(time_SO_Normal_RRA, gastrocnemius_muscle_activation_normal_rra_spring.to_numpy(), label='Spring')\n",
    "axs.plot(time_SO_Strong_RRA, gastrocnemius_muscle_activation_normal_rra.to_numpy(), label='No Spring')\n",
    "axs.set_xlabel('Time (s)')\n",
    "axs.set_ylabel('Activation')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bT4OJvLz8-0J"
   },
   "source": [
    "**Questions**:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Tvb9ALd9oJAv"
   },
   "source": [
    "*  Why did the gastrocnemius activation decrease when the ankle spring was added?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rFr_vdRjoKkw"
   },
   "source": [
    "*  Optional: try increasing or decreasing the spring stiffness. How does this affect the gastrocnemius activation? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rC_GQWCuoL8o"
   },
   "source": [
    "*  Optional: are any other muscles affected by adding the ankle spring?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "i7rhPT2id6gg"
   },
   "source": [
    "## 6.7. Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1XK9UgItWQMF"
   },
   "source": [
    "In this tutorial you have learned how to perform a Static Optimization Analysis, to troubleshoot the different problems associated with it by filtering the motion, using RRA, and tuning the reserve and residual actuators. Also, you have learned to ann a passive element to your model and perform a Static Optimization analysis with it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lLWgd8Noun6Y"
   },
   "source": [
    "## 6.8. Useful Links\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "m8M3dlT5xQcD"
   },
   "source": [
    "> **OpenSim Website:** https://opensim.stanford.edu/\n",
    ">\n",
    "> **OpenSim Python Scripting:** https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53085346/Scripting+in+Python\n",
    ">\n",
    "> **OpenSim API Documentation:** https://simtk.org/api_docs/opensim/api_docs/\n",
    "> \n",
    "> **OpenSim Creator Website:** https://opensimcreator.com/\n",
    "> \n",
    "> **SimTK Website:** https://simtk.org/projects/opensim\n",
    "> \n",
    "> **Biomechanics of Movement's Course:** https://www.youtube.com/channel/UCDNGy0KKNLQ-ztcL5h2Z6zA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rFR2Rf9tk_9A"
   },
   "source": [
    "## 6.9 Acknowledgments\n",
    "\n",
    "Thanks to [OpenSimColab](https://simtk.org/projects/opencolab) project [[3]](https://doi.org/10.1080/10255842.2022.2104607) for creating the first OpenSim Conda package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dxMjdukAnH2N"
   },
   "source": [
    "## 6.10. References\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WX2A_iaYulV1"
   },
   "source": [
    "> [1].   Delp, S. L., Anderson, F. C., Arnold, A. S., Loan, P., Habib, A., John, C. T., Guendelman, E., & Thelen, D. G. (2007). **OpenSim: open-source software to create and analyze dynamic simulations of movement.** *IEEE Transactions on Bio-Medical Engineering*, 54(11), 1940–1950. https://doi.org/10.1109/TBME.2007.901024\n",
    "> \n",
    "> [2] Hicks, J. L., Uchida, T. K., Seth, A., Rajagopal, A., & Delp, S. L. (2015). **Is my model good enough? Best practices for verification and validation of musculoskeletal models and simulations of movement.** *Journal of Biomechanical Engineering*, 137(2), 020905. https://doi.org/10.1115/1.4029304\n",
    ">\n",
    "> [3] Mokhtarzadeh, H., Jiang, F., Zhao, S., & Malekipour, F. (2022). **OpenColab project: OpenSim in Google colaboratory to explore biomechanics on the web.** *Computer Methods in Biomechanics and Biomedical Engineering*, 1–9. https://doi.org/10.1080/10255842.2022.2104607"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "opensim",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
